Observation of shock waves in a strongly interacting Fermi gas 
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We study collisions between two strongly interacting atomic Fermi gas clouds. We observe exotic 
nonlinear hydrodynamic behavior, distinguished by the formation of a very sharp and stable density 
peak as the clouds collide and subsequent evolution into a box-like shape. We model the nonlinear 
dynamics of these collisions using quasi-lD hydrodynamic equations. Our simulations of the time- 
dependent density profiles agree very well with the data and provide clear evidence of shock wave 
formation in this universal quantum hydrodynamic system. 



PACS numbers: 03.75.Ss 

Shock waves are of broad interest and can occur in 
nonlinear quantum hydrodynamic systems when regions 
of high density move with a faster local velocity than re- 
gions of low density resulting in large density gradients. 
In the absence of dissipative or dispersive forces large gra- 
dients would eventually lead to a "gradient catastrophe." 
In cold atomic gases, quantum shock waves have been ob- 
served by using slow light methods to create sharp den- 
sity features in a weakly interacting Bose-Einstein con- 
densate (BEG) [ij . Predictions of the density profiles for 
dispersive shock waves propagating in weakly interacting 
BEG's agree with observations, while predictions for dis- 
sipative classical shock waves only reproduce the slowly 
varying features of the density profiles [5]. In BEG's, 
dispersive shock waves produce soliton trains [1,], which 
also have been observed and modeled for rapidly rotating 
BEG's fil and for merging and splitting BEG's 3- While 
quantum shock waves have been studied theoretically in 
BEG's 0,11,0, they are also of recent interest in theories 
of non-equilibrium electron Fermi gases 0, Q • 

In this Letter, we show that strongly interacting atomic 
Fermi gases provide a new universal medium for studies 
of nonlinear hydrodynamics of quantum matter. Strong 
interactions are characterized by a divergent s-wave scat- 
tering length, which is obtained by using a bias mag- 
netic field to tune to a broad coUisional (Feshbach) res- 
onance. In this so-called unitary regime, the chemical 
potential and the pressure are universal functions of the 
density and temperature. In common with a quark-gluon 
plasma a state of matter that existed microseconds 
after the Big Bang, strongly interacting Fermi gases ex- 
hibit features of nearly perfect (low viscosity) hydrody- 
namics, such as anisotropic (elliptic) flow [10|, and serve 
as a paradigm for other exotic quantum hydrodynamic 
Fermi systems, such as nuclear matter 

Our experiments employ a 50:50 mixture of the two 
lowest hyperfine states of ^Li, confined in a cigar-shaped 



GO2 laser trap, and bisected by a blue-detuned beam at 
532 nm, which produces a repulsive potential. The gas is 
then cooled via forced evaporation near a broad Feshbach 
resonance at 834 G [13] ■ After evaporation, the trap 
is adiabatically recompressed to 0.5% of the initial trap 
depth. This procedure produces two spatially separated 
atomic clouds, containing a total of ~ 10^ atoms per 
spin state. In the absence of the blue-detuned beam, 
the trapping potential is cylindrically symmetric with a 
radial trap frequency of lu^ — ojy = oj± = 27r x 437 Hz 
and an axial trap frequency of Uz 
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27r X 27.7 Hz, where the axial frequency of the optical 
trap is ujoz = 27r x 18.7 Hz and lomz — 2tt x 20.4 Hz 
arises from curvature in the bias magnetic field. When 
the repulsive potential is abruptly turned off, the two 
clouds accelerate toward each other and collide in the 
GO2 laser trap. After a chosen hold time, the GO2 laser 
is turned off, allowing the atomic cloud to expand for 1.5 
ms, after which it is destructively imaged with a 5 /iS 
pulse of resonant light. 

Fig. [T] shows false color absorption images for a colli- 
sion of the atomic clouds at different times after the blue- 
detuned beam is extinguished. Two distinctive features 
are clearly seen in this data: (i) the formation of a central 
peak, which is well-pronounced and robust; (ii) the evo- 
lution of this peak into a box-like shape with very sharp 
boundaries. The observed large density gradients provide 
strong evidence of shock wave formation in this system, 
where the sharp boundaries of the "box" are identified 
as shock wave fronts. Numerical modeling of the hydro- 
dynamic theory for one dimensional motion is used to 
predict the evolution of the atomic density, yielding pro- 
files in good agreement with the data. 

For simplicity, we assume that the cloud is a unitary 
Fermi gas at zero temperature, i.e., we model the cloud 
as a single fluid, consistent with our measurements of the 
sound velocity jT§|. In this case, the local chemical po- 
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FIG. 1: Collision between two unitary Fermi gas clouds in a cigar-shaped optical trap. The clouds are initially separated by 
a repulsive 532 nm optical beam. After the 532 nm beam is extinguished (0 ms), the clouds approach each other. False color 
absorption images show the spatial profiles versus time. Initially, a sharp rise in density occurs in the center of the collision 
zone. At later times the region of high density evolves from a "peak- like" shape into a "box- like" shape. The well defined edges 
of the central zone in the last three images provide evidence of shock wave formation in the unitary Fermi-gas. 



tential has the universal form ^{n^o) = (1 + P)^F{n3D), 
where tpin^D) = |^(37''^'T-3d)^^'^ is the ideal gas lo- 
cal Fermi energy corresponding to the three-dimensional 
density n^D- Here, (3 = —0.61 is a universal scale fac- 
tor [iaIiiJii- 

Neglecting viscous forces, the dynamics for the density 
^3D(r, i) and the velocity field v(r,i), are described by 
the continuity equation. 



dtnsD + V • {nzD^) = 
and the Euler equation. 
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where we assume irrotational flow. Here Utrapi^) — 
imw^r^ -|- imw^z^ is the confining harmonic potential 
of the cigar-shaped trap. 

To determine the initial density profile for the sepa- 
rated clouds, we consider the equilibrium 3D density of 
the Fermi gas in the trap, including a knife-shaped re- 
pulsive blue-detuned beam potential Vrep{z)- The blue- 
detuned laser beam is shaped by a cylindrical lens tele- 
scope, so that the spot size is small compared to the long 
dimension of the cigar-shaped cloud and large compared 
to the transverse dimension. Therefore, the repulsive po- 
tential can be considered to vary only in the z (axial) 
direction according to Vrep{z) — Voexp {—{z — z^Y /tjl). 
We measure the width cr^ = 21.2 /im. The offset zo = 
5 /xm of the focus from the center in the long direction of 
the optical trap is determined by a fit to the first den- 
sity profile at ms. Using the beam intensity and the 
ground state static polarizability of ^Li at 532 nm, we 
find Vb = 12.7 /xK. The initial density profile is then 
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where h = [(2TO/XG/?i^)/(l + /3)]3/V(37r2). In Eq. H 

Rz^i_ — ^2[iQ I (muj]^ _|_) and [iq is the global chemical 

potential, which is determined by normalizing the inte- 
gral of the 3D density to the total number iV of atoms in 
both spin states. For TV = 2 x 10^, we find [Lq = 0.53 /iK, 
i?2 = 220 /im, and i?^ ~ 14 /im. 

We note that yi,Ql[huj\}j — 27, which means that the 
typical number of filled energy levels of transverse quanti- 
zation is large. Therefore, in this paper we use 3D hydro- 
dynamics, Eqs. [l]and[2l and neglect effects of transverse 
quantization even though they are more pronounced in 
regions with lower density. 

We model the dynamics for the one-dimensional mo- 
tion in the long direction of the cigar-shaped trap. Just 
after the blue-detuned beam is extinguished, the initial 
ID density profile is determined by integrating risD of 
Eq. [3]over the transverse dimension r, 



n\i){z) 



27r 



z 

m 



fJ-G 



(4) 



In the following we assume that during the evolution 
the r dependence of Eq. [3] is preserved with the effective 
size of the cloud being a slow function of z and t. We also 
assume that the hydrodynamic velocity is along z axis 
and does not depend on r. Then the subsequent time 
evolution of the density follows the quasi-lD nonlinear 
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hydrodynamic equations: 



dtti — ~dz (nv) 



dfV = -5, — 
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where C 



+ and h 



yjfi/ {mio^) is the oscihator length. For brevity, we have 
omitted the subscript IZ? in Eqs.[5]and|ni The last "vis- 
cosity" term in Eq. [S]is added phenomenologically to de- 
scribe dissipative effects. For the unitary ID fluid, v is 
the effective kinematic viscosity, which has a natural scale 
h/m. It is the only fitting parameter in the theory [l6| . 

For sound wave experiments with a small pulse Vq, one 
can linearize the differential equations ^ and ^ around 
an equilibrium density configuration uq (z) in a harmonic 
trap. Defining n{z,t) = no(z) -I- Sn[z,t), the linearized 
evolution equation for Sn{z,t) (neglecting viscosity) is 
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For a flat background density, i.e., constant uq, with 
= Cuq^^, Eq. [7] reduces to dfSn — c^d'^Sn with the 
sound velocity c = ^2/LtG/5m, in agreement with previ- 
ous theory [l3, [3] and experiment jist . 

To compare the numerical solutions of Eqs. [5] and [6] 
with experiment, we note that the images are taken after 
an additional free expansion for 1.5 ms, during which uid 
continues to slowly evolve in the axial potential of the 
bias magnetic field, i.e., uJz — >■ ujmz = 2tt x 20.4 Hz. We 
assume that during this expansion, the transverse density 
profiles keep the same form, but the radius increases with 
time. Then n3u{r,z) 'n3jy{r/b±,z)/b'^, where b±{t) is 
a transverse scale factor, which obeys b± = tu^bj^^^, 
with &i(0) = 1 and &i(0) [i3,[il[23l- Since the 3D 
pressure scales as rig^^, the ID pressure scales as b_^^^ . 
This leads to a simple modification of Eq. |6l C — )• C(i) = 

c/bT{t). 

We numerically integrate Eqs. [S] and [S] using the mea- 
sured values of the trap frequencies, atom number, and 
the offset, depth, and width of the repulsive potential. 
In the numerical simulation we create and load a den- 
sity array as well as a velocity array with grid spacing 
5z- The initial velocity is set to zero. The simulation 
then updates the density and velocity field in discrete 
time steps 6t according to Eqs. [5] and El The ID den- 
sity profiles are calculated as a function of time after the 
repulsive potential is extinguished. Fig. [2] shows the pre- 
dictions and the data, which are in very good agreement. 
For the simulation curves shown in the figure we use a 
grid of 150 points. To check for numerical consistency, 



we also employ a smoothed-particle- hydrodynamics [21 1 
approach, where the fluid is described by discrete pseudo- 
particles. The results obtained indeed coincide with the 
discretized-grid approach described above. 
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FIG. 2: ID density profiles divided by the total number of 
atoms versus time for two colliding unitary Fermi gas clouds. 
The normalized density is in units of 10~'^//xm per particle. 
Red dots show the measured ID density profiles. Black curves 
show the simulation, which uses the measured trap parame- 
ters and the number of atoms, with the kinetic viscosity as 
the only fitting parameter. 



As shown in Fig. [21 we observe a dramatic evolution 
for the density of the gas. During the collision, a distinct 
and stable density peak forms at the point of collision 
in the center of the trap 



22j . The density gradient at 



the side of the central peak increases from its onset until 
~ 3 ms, at which point the gradient reaches its maximum 
value. A large gradient at the edge of the collision zone 
is maintained throughout the rest of the experiment. For 
most of the data, we find relatively small deviations from 
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the simulation. The largest deviation occurs at 4 ms, 
where the maximum density of the observed central peak 
exceeds that of the simulation by ~ 20%. 

The steep density gradients observed in Fig. [1] suggest 
shock wave formation. A deeper analysis of the simula- 
tion curves provides additional evidence for shock waves. 
Without any dissipation, the numerical integration of the 
quasi- ID theory breaks down due to a "gradient catastro- 
phe." We find that the dissipative force in Eq.[6l which is 
described by the kinematic viscosity coefficient is re- 
quired to attenuate the large density gradients and avoid 
gradient catastrophe. For the data shown in Fig. [2l we 
find that the best fits are obtained with the viscosity pa- 
rameter 1/ = 10 Ti/m. For smaller values of the simula- 
tion produces qualitatively similar results to those shown 
in the figure, only with steeper density gradients at the 
edges of the collision zone. The dissipative term oc v has 
a relatively small effect on the density profiles, unless we 
are in a shock wave regime, where the density gradients 
are large. Hence, the numerical model suggests that the 
large density gradient observed at the edge of the colli- 
sion zone is the leading edge of a dissipative shock wave. 

Our data for a strongly interacting Fermi gas are very 
well described by a simple one-dimensional model based 
on dissipative nonlinear quantum hydrodynamics. In this 
sense, a unitary Fermi gas is different than a weakly inter- 
acting Bose-Einstein condensate, where dispersive terms 
play a major role. 

We derive the one dimensional model, with an effective 

2/5 

one-dimensional chemical potential ~ C n^'j^ , as- 
suming a single fluid near the ground state. However, we 
expect that at higher temperatures, even in the normal 
fluid regime, rapid coUisional equilibrium in the unitary 
gas will produce nearly adiabatic evolution with a three- 
dimensional pressure oc n^/^, and hence an identical 
power-law dependence for the effective one-dimensional 
chemical potential. 

Our observations suggest that collisions of two clouds 
of ultra-cold atoms in the unitary regime are accompa- 
nied by shock wave formation. In future work, it will 
be interesting to study the origin of the effective viscos- 
ity and the effects of transverse quantization. Further, 
the large density gradients produced in the experiments 
suggest that it may be possible to investigate the effects 
of higher derivative dispersive terms in the stress ten- 
sor [2, 0, B [S^l- Finally, the radial density variations 
observed in the two dimensional image is not captured 
in the one-dimensional profiles, but may be studied by 
expanding the analysis to higher dimensions. 
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